Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny

Sara Acevedo

Plan

El libro se puede encontrar acá

Son 16 capítulos
El club solo contempla los 10 primeros
Horario: lunes de 15 a 16 cada dos semanas (1 o 2 capítulos por sesión)

Capítulo 1: Geospatial health

Geospatial health data
Disease mapping
Communication of results

Capítulo 1: Geospatial health data

Los datos en salud permiten detectar problemas de salud pública
Mejoran la eficacia en la respuesa y ayudan a prevenir y controlar enfermedades
El análisis de datos en salud implica distintos métodos, centrados en patrones y conclusiones estadísticas
El análisis espacial usando herramientas SIG juega un papel clave

Capítulo 1: Disease mapping

El mapeo de enfermedades proporciona una visualización de patrones espaciales de las enfermedades
Los mapas resumen zonas de alto riesgo, permiten formular hipótesis y asignar recursos de forma eficaz
Usando modelos jerárquicos bayesianos se puede estimar el riesgo

Capítulo 1: Disease mapping

Métodos bayesianos mencionados son dos:

Markov chain Monte Carlo (MCMC)
Integrated nested Laplace approximation (INLA)

Capítulo 1: Disease mapping

Estos métodos consideran variables y correlación espacial, expresando incertidumbre en la estimación de riesgo
La agregación de datos permite confidencialidad (por ejemplo comunas)
Datos y covariables de alta resolución permiten realizar estimaciones de riesgo precisas

Capítulo 1: Communication of results

La difusión de la información en temas de salud es clave
R ofrece herramientas de comunicación eficaces, como mapas interactivos (leaflet) y gráficos de series temporales (dygraphs)
Permite crear informes reproducibles (RMarkdown ahora Quarto), web interactivas (flexdashboard) y aplicaciones web (Shiny)
La interpretación de estos resultados ayuda a asignar recursos de forma eficiente.

Capítulo 2: Spatial data and R packages for mapping

Types of spatial data
Coordinate reference systems
Shapefiles
Making maps with R

Capítulo 2: Types of spatial data

Tres características: atributo, ubicación de la observación y el dominio

Los tipos de datos según el dominio que representan.

Areal data 🗾
Geostatistical data 🌎 📊
Point patterns 📍

Capítulo 2: Areal data 🗾

El dominio es fijo (de forma regular o irregular) y se divide en un número finito de unidades de área con límites bien definidos.

Capítulo 2: Areal data 🗾

Sudden infant deaths in North Carolina in 1974.

Capítulo 2: Geostatistical data 🌎 📊

En los datos geoestadísticos, el dominio es un conjunto fijo continuo. Es importante señalar que la continuidad se refiere al dominio, y el atributo puede ser continuo o discreto.

Capítulo 2: Geostatistical data

Average rainfall measured at 143 recording stations in Paraná state, Brazil.

Capítulo 2: Point patterns 📍

A diferencia de los datos geoestadísticos y de área, el dominio es al azar.

Capítulo 2: Point patterns 📍

John Snow’s map of the 1854 London cholera outbreak

Capítulo 2: Coordinate reference systems 🌐

Se debe conocer en qué proyección espacial se encuentran los datos.
Las proyecciones están determinadas por un CRS (coordinate reference system).
Cada CRS puede estar definido por un EPSG (sigla de European Petroleum Survey Group), se pueden ver acá
Si quiero hacer operaciones espaciales, todas las capas deben tener el mismo sistema de proyección.

Capítulo 2: Coordinate reference systems 🌐

En Chile, los sistemas de proyección más usados son:

EPSG:4326: WGS 84
EPSG:32719: WGS 84 / UTM zone 19S
EPSG:3857: WGS84 Web (Pseudo)Mercator (Auxiliary Sphere)
EPSG:24879: PSAD56 / UTM zone 19S

Estas slides se basaron en la presentación de Stephanie Orellana para Rladies Madrid

Ejemplo libro (modificado)

Se crea un dataframe al azar con long y lat

d <- data.frame(long = rnorm(3, 0, 1), lat = rnorm(3, 0, 1))
d

       long        lat
1 -1.651048 -0.5910935
2 -2.560355  0.8331618
3  1.302916 -0.5013650

Define las coordenadas

coordinates(d) <- c("long", "lat")
d

class       : SpatialPoints 
features    : 3 
extent      : -2.560355, 1.302916, -0.5910935, 0.8331618  (xmin, xmax, ymin, ymax)
crs         : NA

Ejemplo libro (modificado)

Asigna una proyección (CRS WGS84)

proj4string(d) <- CRS("+proj=longlat +ellps=WGS84
                      +datum=WGS84 +no_defs")
d

class       : SpatialPoints 
features    : 3 
extent      : -2.560355, 1.302916, -0.5910935, 0.8331618  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs

Ejemplo libro (modificado)

Reprojecta desde long/lat a UTM zone 35 south

d_new <- spTransform(d, CRS("+proj=utm +zone=35 +ellps=WGS84
                      +datum=WGS84 +units=m +no_defs +south"))
d_new

class       : SpatialPoints 
features    : 3 
extent      : -2946545, -2461046, 9925476, 10105982  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs

Ejemplo libro (modificado)

Crea nuevas columnas llamadas UTMx y UTMy

d_new$UTMx <- coordinates(d_new)[, 1]
d_new$UTMy <- coordinates(d_new)[, 2]

Capítulo 2: Shapefiles 🗺️

Shapefiles almacenan ubicación, forma y atributos de puntos, líneas y polígonos
Un shapefile no es un archivo único, sino que consiste en una colección de archivos relacionados que tienen diferentes extensiones con un nombre común
Un shapefile tiene tres archivos obligatorios con extensiones .shp, .shx y .dbf

Capítulo 2: Shapefiles 🗺️

En R, podemos leer shapefiles utilizando la función readOGR() del paquete rgdal, o también la función st_read() del paquete sf

Danger

during October 2023 rgdal, rgeos and maptools will be archived on CRAN, and packages with strong dependencies on the retiring packages must be either upgraded to use sf, terra or other alternatives or work-arounds by or before that time.

Capítulo 2: Shapefiles 🗺️

nameshp <- system.file("shape/nc.shp", package = "sf")

map <- st_read(nameshp, quiet = TRUE)

class(map)

[1] "sf"         "data.frame"

Capítulo 2: Shapefiles 🗺️

plot(map)

Capítulo 2: Making maps with R 🗺️

Los paquetes mas útiles para hacer mapas son

ggplot2
leaflet
mapview
tmap

Capítulo 2: Making maps with R 🗺️

Acá usaré el mismo código del libro pero con datos de CEDEUS. Primero cargo los datos

uso_suelo<- read_sf("https://opendata.arcgis.com/datasets/08dea558cea94b64bee9074ed0fbfd4f_0.geojson")

class(uso_suelo) # ya es sf, sino se transforma con st_as_sf()

[1] "sf"         "tbl_df"     "tbl"        "data.frame"

Capítulo 2: Making maps with R 🗺️ ggplot2

ggplot(uso_suelo) + 
  geom_sf(aes(fill = uso_suelo)) +
  theme_bw()

Capítulo 2: Making maps with R 🗺️ ggplot2

Puedo usar escalas de colores definidas

ggplot(uso_suelo) + 
  geom_sf(aes(fill = Shape__Area)) +
  scale_fill_viridis() +
  theme_bw()

Capítulo 2: Making maps with R 🗺️ leaflet

Leaflet es una libreria de JavaScript de código abierto para mapas interactivos
Leaflet requiere una proyección específica: EPSG code 4326

Capítulo 2: Making maps with R 🗺️ leaflet

Revisamos la proyección

st_crs(uso_suelo)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Capítulo 2: Making maps with R 🗺️ leaflet

Transformamos a EPSG code 4326

uso_suelo_leaflet<- st_transform(uso_suelo, 4326)

st_crs(uso_suelo_leaflet)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

Capítulo 2: Making maps with R 🗺️ leaflet

pal <- colorNumeric("BrBG", domain = uso_suelo$Shape__Area)
leaflet(uso_suelo_leaflet) %>%
  addTiles() %>%
  addPolygons(
    color = "white", fillColor = ~ pal(Shape__Area),
    fillOpacity = 1
  ) %>%
  addLegend(pal = pal, values = ~Shape__Area, opacity = 1)

Capítulo 2: Making maps with R 🗺️ mapview()

mapview(uso_suelo, zcol = "uso_suelo")

Capítulo 2: Making maps with R 🗺️ mapview()

Puedo usar escalas de colores de otros paquetes

library(RColorBrewer)
mapview(uso_suelo,
  zcol = "uso_suelo",
  map.types = "CartoDB.DarkMatter",
  col.regions = colorRampPalette(brewer.pal(7, "Paired")) 
)

Capítulo 2: Making maps with R 🗺️ mapview()

Tiene muchas opciones, por ejemplo mapas al lado de otro

uso <- mapview(uso_suelo, zcol = "uso_suelo")
area <- mapview(uso_suelo, zcol = "Shape__Area")
mapa <- leafsync::sync(uso, area)
mapa

Capítulo 2: Making maps with R 🗺️ tmap

tm_shape(uso_suelo) +
  tm_polygons("uso_suelo")

Conclusiones

Datos de salud 🤝 Datos espaciales
Antes de usar los datos: revisar proyección 🌐
Revisar paquetes y su documentación (rgdal y rgeos 🥶)